Survey Data Analysis Case Study

# get Pronoun * PSA interaction from Exp2 production model
load("r_data/exp2.RData")

exp2_r_effect_size <- exp2_m_prod@model |>
  tidy() |>
  filter(term == "Pronoun=They_HeShe:PSA=GenLang") |>
  pull(estimate) |>
  round(2)

exp2_r_effect_size       # log-odds
exp(exp2_r_effect_size)  # odds ratio

# start with 108 participants each doing 30 trials
exp3_pw_data_struct <- data.frame(
  Participant = rep(as.factor(1:108), each = 30),
  Trial = rep(as.factor(1:30), 108)
)

# Trials are split between 3 Pronoun Pair conditions, which are contrast-coded
# to compare:
# (1) They|HeShe vs HeShe|They + HeShe|SheHe
# (2) HeShe|They vs HeShe\|SheHe
exp3_pw_data_struct <- exp3_pw_data_struct |>
  bind_cols(
    "Pronoun" = rep(rep(factor(c("He", "She", "They")), each = 10), 108)
  )
contrasts(exp3_pw_data_struct$Pronoun) <- cbind(
  "_T vs HS" = c(.33, .33, -.66),
  "_H vs S"  = c(-.5, .5, 0)
)

# Nametag and Introduction conditions vary in a 2x2 between-P design, and both
# are mean-centered effects coded.
exp3_pw_data_struct <- exp3_pw_data_struct |>
  bind_cols(
    "Nametag" = rep(rep(factor(c(0, 0, 1, 1)), each = 30), 108 / 4),
    "Intro" = rep(rep(factor(c(0, 1, 0, 1)), each = 30), 108 / 4)
  )

contrasts(exp3_pw_data_struct$Nametag) <- cbind("_No_Yes" = c(-.5, .5))
contrasts(exp3_pw_data_struct$Intro) <- cbind("_No_Yes" = c(-.5, .5))

# Item is defined as each unique image-name-pronoun combination. There are 6
# sets of characters, and each list sees 3, making 18 unique characters.
exp3_pw_data_struct <- exp3_pw_data_struct |>
  bind_cols(
    "Character" = rep(as.factor(1:18), each = 30 / 3, 108 / 6)
  )
str(exp3_pw_data_struct)

exp3_pw_data_struct |>
  group_by(Nametag, Intro) |>
  summarise(n_distinct(Participant))

# The closest thing to existing data is the Exp2 (written) production task.
# Since interpreting effect sizes is apparently more complicated for logistic
# regression, let's go with the Exp2 results as a baseline. That's a rough
# estimate of how much harder they/them is to produce than he/him and she/her.
# And let's set the hypothetical Nametag and Introduction effects to be about
# the same size as the PSA. Hopefully that's small enough to be kind of
# conservative with the power analysis, but not aiming for effects too small to
# be practically relevant.
exp2_m_prod_fixed <- exp2_m_prod@model |>
  tidy() |>
  filter(effect == "fixed") |>
  select(term, estimate)
exp2_m_prod_fixed

# Predictions for Exp3 based on ranges from Exp2:
exp3_pw_fixed <- c(
  +0.75,  # Intercept                    Medium
  +3.00,  # Pronoun: T vs HS             Largest
  -0.10,  # Pronoun: H vs S              NS, maybe small
  +0.10,  # Nametag                      NS, maybe small
  +0.10,  # Introduction                 NS, maybe small
  -2.00,  # Pronoun: T vs HS * Nametag   Same size as PSA interaction
  -0.10,  # Pronoun: H vs S  * Nametag   NS, maybe small
  -2.00,  # Pronoun: T vs HS * Intro     Same size as PSA interaction
  -0.10,  # Pronoun: H vs S  * Intro     NS, maybe small
  +0.25,  # Nametag * Intro              Maybe small
  -2.00,  # 3 way T vs HS                Same size as PSA interaction
  -0.10   # 3 way H vs S                 NS, maybe small
)

# The model for the Exp2 production task only converged with random intercepts
# by item, and no random effects by participant.
exp2_m_prod_random <- VarCorr(exp2_m_prod@model)

# The model for the Exp1 production task only converged with random intercepts
# and slopes by participant, and no random effects by item.
load("r_data/exp1.RData")
exp1_m_prod_random <- VarCorr(exp1a_m_prod@model)


# So, I'll combine those two as a starting place to estimate the random effects.
# It's possible the actual data won't converge with the maximal random effects
# structure, but for now let's assume it will.
exp3_pw_random <- exp1_m_prod_random
exp3_pw_random[["Item"]] <- exp2_m_prod_random[["Name"]]

# Create model with this data structure, fixed effects, and random effects
exp3_pw_m_108 <- makeGlmer(
  formula = SimAcc ~ Pronoun * Nametag * Intro +
    (Pronoun | Participant) + (1 | Character),
  family = binomial,
  fixef = exp3_pw_fixed,
  VarCorr = exp3_pw_random,
  data = exp3_pw_data_struct
)
summary(exp3_pw_m_108)

# Simulate data
exp3_pw_sim_data <- doSim(exp3_pw_m_108)
exp3_pw_data_struct <- exp3_pw_data_struct |>
  bind_cols("SimAcc" = exp3_pw_sim_data)

summary(exp3_pw_data_struct)

# Code to run simulation:
powerSim(
  exp3_pw_m_108,
  nsim = 1000,
  test = fixed("Pronoun_T vs HS:Nametag_No_Yes", "z")
)

# Then extend model to larger N
exp3_pw_m_132 <- extend(exp3_pw_m_108, along = "Participant", n = 132)

# Load and join results
exp3_pw_results <- bind_rows(
    .id = "sim",
    "2_108" = readRDS("r_data/exp3_power_2way_N108.RDA") |> summary(),
    "2_132" = readRDS("r_data/exp3_power_2way_N132.RDA") |> summary(),
    "2_156" = readRDS("r_data/exp3_power_2way_N156.RDA") |> summary(),
    "2_180" = readRDS("r_data/exp3_power_2way_N180.RDA") |> summary(),
    "3_132" = readRDS("r_data/exp3_power_3way_N132.RDA") |> summary(),
    "3_156" = readRDS("r_data/exp3_power_3way_N156.RDA") |> summary()
  ) |>
  mutate(
    n_participants = str_sub(sim, 3),
    effect = case_when(
      str_sub(sim, 0, 1) == "2" ~ "Pronoun * Nametag/Intro",
      str_sub(sim, 0, 1) == "3" ~ "Pronoun * Nametag * Intro"
    )
  ) |>
  column_to_rownames(var = "sim")

# Options:
#   s   sync data from AWS
#   p   check participant list
#   z   unzip and sort audio files
#   c   run tests on PCIbex output and audio file names
#   t   transcribe


while getopts "spzct" option; do
  case $option in
    s)  # Sync audio data from S3
        echo "Getting data from AWS"
        cd ../data/s3/
        aws s3 sync s3://they3 .
        cd ../../preprocessing/
        ;;
    p)  # Get list of participants from PCIbex data to update participant list
        echo "Checking audio data to see what needs to be added to the participant list"
        Rscript participant_list.R
        ;;
    z)  # Unzip the audio data and convert it to WAV files in dirs for each participant
        echo "Unzipping, converting, and sorting the audio files"
        python s3_to_wav.py
        ;;
    c)  # Check output
        echo "Checking the audio file names against the PCIbex data"
        Rscript check_output.R
        ;;
    t)  # Transcribe
        echo "Transcribing"
        python transcribe.py
        ;;
  esac
done

import os
import whisper
import pandas as pd
from pathlib import Path


# ---- Helper functions ----- #
def make_transcription_df():
    """Set up df for transcription data.

    Returns:
        df: columns for `participant_id`, `prolific_id`, `trial_id`, indexed by
            `file_path`
    """    
    transcriptions = []
    participant_dirs = [
        p for p in audio_dir.iterdir()
        if not p.match("*temp*") and not p.match("*incomplete*")
    ]

    for p in participant_dirs:
        audio_list = [a.stem for a in p.glob('*.wav')]
        trials = [get_trial_info(p.name, a) for a in audio_list]
        df = pd.DataFrame(
            trials,
            columns=['file_path', 'participant_id', 'prolific_id', 'trial_id'],
        )
        df = df.set_index('file_path').sort_values(by='trial_id')
        transcriptions.append(df)
    return transcriptions


def get_trial_info(p_dir, file_name):
    """Get trial info from the name of the audio file.

    Args:
        p_dir (str): dir for participant data
        file_name (str): audio file within participant's data dir 

    Returns:
        list: .wav file name (Path), participant ID (str), prolific ID (str),
            and trial ID (str)
    """    
    participant_id, prolific_id = p_dir.split('_')
    trial_id = file_name.removeprefix(prolific_id + '_').removesuffix('.wav')
    return [
        audio_dir / p_dir / f"{file_name}.wav",
        participant_id, prolific_id, trial_id
    ]


def run_whisper_on_participant(df, model):
    """Use whisper to transcribe a trial.

    Args:
        df (df): structure for transcription data from `make_transcription_df()`,
            which has `participant_id` as the first column and is indexed by the
            path to the audio file
        model (whisper model): whisper model loaded (using small English-only)
    """    
    participant_id = df.iloc[0, 0]
    file_path = text_dir / f"{participant_id}_whisper.csv"
    if not os.path.exists(file_path):
        print(participant_id)
        df['text'] = df.index.map(lambda t: whisper.transcribe(model, str(t))['text'])
        print(df['text'])
        df.to_csv(os.path.join(file_path))
     

# ---- Main function ----- #
def transcribe_trials():
    """Main function to transcribe .wav files using whisper."""    
    transcriptions = make_transcription_df()
    model = whisper.load_model('medium.en')
    for p in transcriptions:
        run_whisper_on_participant(p, model)
        
    return transcriptions


# ---- Run ----- #
audio_dir = Path('..') / 'data' / 'exp2_audio'
text_dir = Path('..') / 'data' / 'exp2_transcription'

transcribe_trials()

library(here)
library(tidyverse)
library(testthat)
library(readxl)


# Load data from coded CSVs----
df <- list.files(
  path = "data/exp2_coding",
  pattern = "*coded.csv",
  full.names = TRUE
) |>
  set_names() |>
  map(read.csv) |>
  list_rbind(names_to = "participant_id") |>
  mutate(participant_id = str_sub(str_split_i(participant_id, "/", 3), 0, 4)) |>
  mutate(
    .before = trial_id,
    character_list = str_remove(str_split_i(trial_id, "_", 2), "list")
  ) |>
  mutate(across(
    c(ends_with("id"), contains("pronoun"), condition, character_list),
    as.factor
  )) |>
  mutate(across(
    starts_with("transcription"),
    \(x) {
      case_when(
        x == "" ~ NA,
        x == " ." ~ NA,
        .default = str_replace(str_trim(x), "Jamie", "Jaime")  # fix spelling
      )
    }
  )) |>
  mutate(drop_trial = case_when(  # don't drop trial if sentence is cut off
    str_detect(transcription_manual, "Alex") ~ 0,
    str_detect(transcription_manual, "Casey") ~ 0,
    str_detect(transcription_manual, "Jaime") ~ 0,
    str_detect(transcription_manual, "Jordan") ~ 0,
    str_detect(transcription_manual, "Sam") ~ 0,
    str_detect(transcription_manual, "Taylor") ~ 0,
    .default = drop_trial
  )) |>
  filter(str_starts(participant_id, "P"))  # drop incompletes


# Calculate pronoun produced and accuracy----
df2 <- df |>
  mutate(
    sum = his + her + their,
    multiple_pronouns = ifelse(sum > 1, 1, 0),
    he_loc = ifelse(
      multiple_pronouns == 1 & he == 1,
      str_locate(str_to_lower(transcription_manual), "\\bhe"), 0
    ),
    his_loc = ifelse(
      multiple_pronouns == 1 & his == 1,
      str_locate(str_to_lower(transcription_manual), "\\bhi"), 0
    ),
    she_loc = ifelse(
      multiple_pronouns == 1 & she == 1,
      str_locate(str_to_lower(transcription_manual), "\\bshe\\b"), 0
    ),
    her_loc = ifelse(
      multiple_pronouns == 1 & her == 1,
      str_locate(str_to_lower(transcription_manual), "\\bher\\b"), 0
    ),
    they_loc = ifelse(
      multiple_pronouns == 1 & they == 1,
      str_locate(str_to_lower(transcription_manual), "\\bthey\\b"),
      0
    ),
    their_loc = ifelse(
      multiple_pronouns == 1 & their == 1,
      str_locate(str_to_lower(transcription_manual), "\\bthei"),
      0
    ),
    max_loc = ifelse(
      multiple_pronouns == 1,
      pmax(he_loc, his_loc, she_loc, her_loc, they_loc, their_loc),
      NA
    ),
    pronoun_produced = factor(
      case_when(
        sum == 0 ~ "none",
        multiple_pronouns == 0 & (his == 1 | he == 1) ~ "his",
        multiple_pronouns == 0 & (her == 1 | she == 1) ~ "her",
        multiple_pronouns == 0 & (their == 1 | they == 1) ~ "their",
        multiple_pronouns == 1 & max_loc == he_loc ~ "his",
        multiple_pronouns == 1 & max_loc == his_loc ~ "his",
        multiple_pronouns == 1 & max_loc == she_loc ~ "her",
        multiple_pronouns == 1 & max_loc == her_loc ~ "her",
        multiple_pronouns == 1 & max_loc == they_loc ~ "their",
        multiple_pronouns == 1 & max_loc == their_loc ~ "their",
        is.na(max_loc) & participant_id == "P329" &
          str_detect(transcription_manual, "chocolate") ~ "his",
        is.na(max_loc) & participant_id == "P492" &
          str_detect(transcription_manual, "pencil") ~ "their"
      ),
      levels = c("his", "her", "their", "none")
    ),
    accuracy = case_when(
      pronoun_produced == "none" ~ NA,
      as.character(pronoun_produced) == "his" &
        as.character(target_pronoun) == "he" ~ 1,
      as.character(pronoun_produced) == "her" &
        as.character(target_pronoun) == "she" ~ 1,
      as.character(pronoun_produced) == "their" &
        as.character(target_pronoun) == "they" ~ 1,
      .default = 0
    )
  )

# Drop trials with no data, drop calculation/extra columns, add item ID----
df3 <- df2 |>
  filter(drop_trial == 0 & !is.na(transcription_manual)) |>
  rename(transcription = transcription_manual) |>
  select(
    -prolific_id, -correct_description, -transcription_whisper, -name_only,
    -drop_trial, -notes, -sum, -ends_with("loc")
  )

character_lists <- read_excel(
  path = here("materials", "exp2_stimuli.xlsx"),
  sheet = "Character Sets",
  skip = 2,
  col_types = c("skip", rep("text", 7), rep("skip", 9))
) |>
  rename_with(str_to_lower) |>
  mutate(across(everything(), as.factor))

df3 <- df3 |>
  left_join(
    character_lists |> select(list, item, pronouns),
    by = join_by(character_list == list, target_pronoun == pronouns),
    relationship = "many-to-one"
  ) |>
  relocate(item, .after = target_pronoun) |>
  rename(target_id = item) |>
  select(-character_list)


# Checks----
test_that("No NAs in pronoun variables", {
  expect_false(any(is.na(df3$he)))
  expect_false(any(is.na(df3$his)))
  expect_false(any(is.na(df3$she)))
  expect_false(any(is.na(df3$her)))
  expect_false(any(is.na(df3$they)))
  expect_false(any(is.na(df3$their)))
  expect_false(any(is.na(df3$name_only)))
})

test_that("Pronoun variables match regex", {
  df_test <- df3 |>
    mutate(transcription = str_to_lower(transcription)) |>
    select(participant_id, transcription, he, his, she, her, they, their) |>
    mutate(
      has_he = case_when(
        str_detect(transcription, "\\bhe\\b") ~ 1,
        .default = 0
      ),
      has_his = case_when(
        str_detect(transcription, "\\bhis\\b") ~ 1,
        str_detect(transcription, "\\bhi-") ~ 1,
        .default = 0
      ),
      has_she = case_when(
        str_detect(transcription, "\\bshe\\b") ~ 1,
        .default = 0
      ),
      has_her = case_when(
        str_detect(transcription, "gave her a") ~ 0,
        str_detect(transcription, "her glasses back") ~ 0,
        str_detect(transcription, "pencil to h--") ~ 1,
        str_detect(transcription, "\\bher\\b") ~ 1,
        .default = 0
      ),
      has_they = case_when(
        str_detect(transcription, "\\bthey\\b") ~ 1,
        .default = 0
      ),
      has_their = case_when(
        str_detect(transcription, "\\btheir\\b") ~ 1,
        str_detect(transcription, "to thei-") ~ 1,
        str_detect(transcription, "to theirs") ~ 1,
        str_detect(transcription, "chocolate to th--") ~ 1,
        .default = 0
      )
    )

  expect_true(all(df_test$he == df_test$has_he))
  expect_true(all(df_test$his == df_test$has_his))
  expect_true(all(df_test$she == df_test$has_she))
  expect_true(all(df_test$her == df_test$has_her))
  expect_true(all(df_test$they == df_test$has_they))
  expect_true(all(df_test$their == df_test$has_their))
})

test_that("All multiple pronoun trials coded as disfluencies and final coded", {
  df_dis <- df3 |>
    filter(multiple_pronouns == 1 & disfluency == 0) |>
    filter(!(he == 1 & his == 1)) |>
    filter(!(she == 1 & her == 1)) |>
    filter(!(they == 1 & their == 1))
  expect_equal(nrow(df_dis), 0)

  expect_false(any(is.na(df3$pronoun_produced)))
})


# Counts----
participants_missing_trials <- df |>
  filter(!is.na(transcription_manual) & drop_trial == 0) |>
  summarise(.by = participant_id, n = n_distinct(trial_id))
sum(participants_missing_trials$n)

df3 |>
  summarise(
    .by = c(condition, target_pronoun),
    mean = mean(accuracy, na.rm = TRUE)
  ) |>
  arrange(target_pronoun)


# Export----
write_csv(df3, file = "data/exp2_pronouns.csv")

library(tidyverse)
library(janitor)
library(readxl)

# MAIN DF----
## Read PCIbex output----
d_survey <- list.files(path = "data/exp2_PCIbex/", full.names = TRUE) |>
  map_df(
    ~read.csv(., header = FALSE, fill = TRUE, col.names = paste("V", 1:26))
  ) |>
  filter(str_detect(V.1, "#") == FALSE) |>  # drop PennController comments
  select(V.13, V.6, V.9, V.15, V.10, V.11) |>  # drop PennController extra cols
  rename(  # name PennController output columns
    trial_type  = V.6,  # trial label in PCIbex
    trial_part  = V.9,  # type of trial data
    parameter   = V.10,
    response    = V.11,
    prolific_id = V.13,  # first trial variable saved is prolific_id
    trial_item  = V.15
  ) |>
  filter(str_detect(  # get demographics and familiarity questions
    trial_type, "demographics|sentences|they|transphobia"
  )) |>
  # Remove status update rows for them
  filter(parameter != "_Header_" & parameter != "_Trial_") |>
  filter(parameter != "First" & parameter != "Unselect") |>
  filter(parameter != "Status" & parameter != "Filename") |>
  filter(!(
    # remove rows that indicate last item selected in check box
    parameter == "Choice" &
      (
        trial_part == "enter_they" | trial_part == "enter_trans" |
          trial_part == "enter_sexuality" | trial_part == "enter_race"
      )
  )) |>
  filter(!(parameter == "Final" & response == ""))  # write-in box empty


## Match to Participant ID----
participant_list <- "data/participant_list.xlsx" |>
  read_xlsx(sheet = 1, range = cell_cols(1:10)) |>
  clean_names() |>
  select(ends_with("id"), condition) |>
  filter(!is.na(participant_id)) |>
  mutate(across(everything(), as.factor))

d_survey <- d_survey |>
  left_join(participant_list, by = "prolific_id") |>
  relocate(participant_id, .before = 1) |>
  select(-prolific_id, -condition)

## Exclusions----
d_survey <- d_survey |> filter(str_detect(participant_id, "P")) |> droplevels()


## Question categories & items----
d_survey <- d_survey |>
  mutate(
    .after = participant_id,
    category =
      case_when(
        trial_type == "rate_sentences" ~ "Sentence Naturalness Ratings",
        trial_type == "transphobia_scale" ~ "Transphobia Scale",
        trial_part == "enter_they" ~ "Familiarity With They/Them Pronouns",
        str_detect(trial_part, "intro|nametag") ~
          "Familiarity With Pronoun-Sharing Practices",
        str_detect(trial_part, "age") ~ "Age",
        str_detect(trial_part, "gender") ~ "Gender",
        str_detect(trial_part, "trans") ~ "Transgender & Gender-Diverse",
        str_detect(trial_part, "sexuality") ~ "Sexuality",
        str_detect(trial_part, "race") ~ "Race/Ethnicity",
        str_detect(trial_part, "english") ~ "English Experience",
        str_detect(trial_part, "ed") ~ "Education"
      ) |>
      as.factor(),
    item =
      case_when(
        trial_part == "enter_intro_others" ~ "Intros: Others",
        trial_part == "enter_intro_self" ~ "Intros: Self",
        trial_part == "enter_nametags_others" ~ "Nametags: Others",
        trial_part == "enter_nametags_self" ~ "Nametags: Self",
        str_detect(parameter, "for myself") ~ "Myself",
        str_detect(parameter, "am close to") ~ "Close To",
        str_detect(parameter, "have met") ~ "Have Met",
        str_detect(parameter, "have not met")  ~ "Heard About",
        str_detect(parameter, "had not heard") ~ "Not Heard About",
        trial_item != "" ~ trial_item,
        str_detect(category, "Trans|Sexuality|Race") ~ parameter,
        str_detect(category, "Ed|Eng") ~ response,
        str_detect(category, "Gender") ~ category
      ) |>
      recode_factor(
        "generic" = "Generic",
        "each" = "Each",
        "every" = "Every",
        "neu" = "Neutral\nName",
        "fem" = "Fem\nName",
        "masc" = "Masc\nName"
      ) |>
      str_replace_all(c(
        "%2C" = ",", "2 year" = "2-year", "4 year" = "4-year", "term:" = "term"
      )) |>
      as.factor()
  ) |>
  select(-starts_with("trial"), -parameter)


## Response types----
d_survey$response <- d_survey$response |>
  str_replace_all(c("%2C" = ",", "2 year" = "2-year", "4 year" = "4-year"))

d_survey <- d_survey |>
  mutate(
    response_num = case_when(
      !is.na(as.numeric(response)) ~ as.numeric(response),
      is.na(as.numeric(response)) ~ NA
    ),
    response_bool = case_when(
      response == "checked" ~ TRUE,
      response == "unchecked" ~ FALSE,
      response != "checked" & response != "unchecked" ~ NA
    ),
    response_cat = case_when(
      is.na(response_num) & is.na(response_bool) ~ response,
      .default = NA
    ),
    item = case_when(
      !is.na(item) ~ item,
      category == "Age" & response_num <= 24 ~ "18–24",
      category == "Age" & response_num >= 25 & response_num <= 34 ~ "25–34",
      category == "Age" & response_num >= 35 & response_num <= 44 ~ "35–44",
      category == "Age" & response_num >= 45 & response_num <= 54 ~ "45–54",
      category == "Age" & response_num >= 55 & response_num <= 64 ~ "55–64",
      category == "Age" & response_num >= 65 & response_num <= 74 ~ "65–74",
      response_num >= 75 ~ "75+"
    )
  ) |>
  mutate(across(where(is.character), as.factor)) |>
  select(-response) |>
  filter(!is.na(item))


## Recode gender----
d_survey |>
  filter(category == "Gender") |>
  pull(response_cat) |>
  droplevels() |>
  unique()

# Group similar responses
d_survey$response_cat <- d_survey$response_cat |>
  recode_factor(
    "female" = "Woman", "f" = "Woman", "Femal" = "Woman", "Female" = "Woman",
    "FEMALE" = "Woman", "woman" = "Woman", "WOMAN" = "Woman", "Female " = "Woman",
    "female/woman" = "Woman", "Female/Woman" = "Woman", "cis woman" = "Woman",
    "Cisfemale" = "Woman", "cisgender woman" = "Woman", "transwoman" = "Woman",
    "male" = "Man", "MALE" = "Man", "Male" = "Man", "Male " = "Man", "Man" = "Man",
    "cis-gender male" = "Man", "cis male" = "Man", "TRANS MAN" = "Man",
    "Transgender Man" = "Man", "Nonbinary" = "Nonbinary spectrum",
    "nonbinary" = "Nonbinary spectrum", "Non-binary" = "Nonbinary spectrum",
    "non binary" = "Nonbinary spectrum", "Transfem nonbinary" = "Nonbinary spectrum",
    "Male and nonbinary" = "Nonbinary spectrum", "she/they" = "Nonbinary spectrum",
    "genderfluid" = "Nonbinary spectrum", "Genderfluid" = "Nonbinary spectrum",
    "questioning" = "Questioning"
  )

d_survey |>
  filter(category == "Gender") |>
  pull(response_cat) |>
  droplevels() |>
  unique()


## Write-in responses----
d_survey$item <- d_survey$item |>
  recode_factor("Final" = "I use a different term")

# Just keep one row for yes to diff term + write-in box
d_survey <- d_survey |>
  mutate(
    response_bool = case_when(
      item == "I use a different term" & !is.na(response_cat) ~ TRUE,
      .default = response_bool
    )
  ) |>
  filter(
    !(item == "I use a different term" & is.na(response_cat) & response_bool == TRUE)
  ) |>
  filter(response_cat != "Normal" | is.na(response_cat))  # asshole response that also checked straight


## Add missing data----
missing <- tibble(
  participant_id = c(
    rep("P277", 31), rep("P278", 31), rep("P419", 31),
    rep("P482", 31), rep("P502", 31)
  ),
  category = rep(
    c(
      "Age", "Education", "English Experience",
      rep("Familiarity With Pronoun-Sharing Practices", 4),
      rep("Familiarity With They/Them Pronouns", 5),
      "Gender", "Race/Ethnicity", "Sexuality",
      rep("Sentence Naturalness Ratings", 6),
      "Transgender & Gender-Diverse",
      rep("Transphobia Scale", 9)
    ),
    5
  ),
  item = rep(
    c(
      rep("Missing Data", 3),
      "Intros: Others", "Intros: Self", "Nametags: Others", "Nametags: Self",
      "Myself", "Close To", "Have Met", "Heard About", "Not Heard About",
      "Missing Data", "Missing Data", "Missing Data",
      "Masc Name", "Fem Name", "Neutral Name", "Generic", "Every", "Each",
      "Missing Data",
      paste(
        "I am uncomfortable around people who don’t conform to traditional",
        "gender roles, e.g., aggressive women or emotional men."
      ),
      "I avoid people on the street whose gender is unclear to me.",
      paste(
        "I think there is something wrong with a person who says that they",
        "are neither a man nor a woman."
      ),
      paste(
        "I would be upset if someone I’d known a long time revealed to me",
        "that they used to be another gender."
      ),
      paste(
        "When I meet someone, it is important for me to be able to identify",
        "them as a man or a woman."
      ),
      "I believe that a person can never change their gender.",
      paste(
        "A person’s genitalia define what gender they are, e.g., a penis",
        "defines a person as being a man, a vagina defines a person as being a",
        "woman."
      ),
      paste(
        "I don’t like it when someone is flirting with me, and I can’t tell",
        "if they are a man or a woman."
      ),
      "I believe that the male/female dichotomy is natural."
    ),
    5
  )
)

missing <- missing |>
  mutate(
    response_num = as.numeric(NA),
    response_bool = NA,
    response_cat = as.character(NA)
  ) |>
  mutate(across(where(is.character), as.factor))

d_survey <- bind_rows(d_survey, missing) |>
  distinct() |>
  mutate(participant_id = as.factor(as.character(participant_id)))


# Aggregates----
## Age----
agg_age <- d_survey |>
  filter(category == "Age") |>
  select(participant_id, response_num) |>
  arrange(participant_id) |>
  rename(age = response_num)


## TGD----
agg_TGD <- d_survey |>
  filter(
    category == "Transgender & Gender-Diverse" & response_bool == TRUE
  ) |>
  select(participant_id, item, response_bool) |>
  mutate(
    response_coded = case_when(
      str_detect(item, "is different") ~ 1,
      item == "I consider myself transgender" ~ 1,
      .default = 0
    )
  ) |>
  summarise(
    .by = participant_id,
    TGD = sum(response_coded) |> recode(`2` = 1)
  )


## LGBQ----
agg_LGBQ <- d_survey |>
  filter(category == "Sexuality" & response_bool == TRUE) |>
  select(participant_id, item, response_bool) |>
  mutate(response_coded = ifelse(str_detect(item, "As|Bi|Gay|Queer"), 1, 0)) |>
  summarise(
    .by = participant_id,
    LGBQ = sum(response_coded) |> recode(`2` = 1, `3` = 1)
  )


## Transphobia scale----
agg_TS <- d_survey |>
  filter(category == "Transphobia Scale") |>
  select(participant_id, response_num) |>
  mutate(response_coded = response_num - 1) |>
  summarise(.by = participant_id, gender_beliefs = sum(response_coded))


## Sentence ratings----
agg_ratings <- d_survey |>
  filter(category == "Sentence Naturalness Ratings") |>
  select(participant_id, item, response_num) |>
  mutate(
    type = ifelse(
      str_detect(item, "Name"), "rating_name", "rating_indefinite"
    )
  ) |>
  summarise(.by = c(participant_id, type), rating = mean(response_num)) |>
  pivot_wider(names_from = "type", values_from = "rating")


## Familiarity with using they/them----
agg_they <- d_survey |>
  filter(str_detect(category, "They/Them") & response_bool == TRUE) |>
  select(participant_id, item, response_bool) |>
  pivot_wider(names_from = item, values_from = response_bool) |>
  mutate(Myself_Close = ifelse(
    Myself == TRUE & `Close To` == TRUE, TRUE, NA
  )) |>
  mutate(
    .keep = c("unused"),
    familiarity = case_when(
      Myself_Close == TRUE | `Close To` == TRUE | Myself == TRUE ~ 3,
      `Have Met` == TRUE ~ 2,
      `Heard About` == TRUE | `Not Heard About` == TRUE ~ 1
    )
  )


## Familiarity with pronoun-sharing----
agg_sharing <- d_survey |>
  filter(str_detect(category, "Sharing")) |>
  select(participant_id, item, response_cat) |>
  mutate(response_coded = case_when(
    response_cat == "Always"    | response_cat == "All"   ~ 5,
    response_cat == "Usually"   | response_cat == "Most"  ~ 4,
    response_cat == "Sometimes" | response_cat == "Some"  ~ 3,
    response_cat == "Rarely"    | response_cat == "A few" ~ 2,
    str_detect(response_cat, "prefer not to") ~ 1,
    str_detect(response_cat, "not heard") ~ 0,
    response_cat == "None"                        ~ 0
  )) |>
  summarise(.by = participant_id, sharing = sum(response_coded))


## Merge---
d_agg <- participant_list |>
  select(participant_id, condition) |>
  left_join(agg_age, by = "participant_id") |>
  left_join(agg_LGBQ, by = "participant_id") |>
  left_join(agg_TGD, by = "participant_id") |>
  left_join(agg_ratings, by = "participant_id") |>
  left_join(agg_sharing, by = "participant_id") |>
  left_join(agg_they, by = "participant_id") |>
  left_join(agg_TS, by = "participant_id")


# Demographics table----
d_demographics <- d_survey |>
  filter(
    category %in% c(
      "Age", "Gender", "Transgender & Gender-Diverse", "Sexuality",
      "Race/Ethnicity", "Education", "English Experience"
    )
  ) |>
  filter(response_bool == TRUE | is.na(response_bool)) |>
  mutate(group = case_when(
    category == "Gender" ~ as.character(response_cat),
    category == "English Experience" ~ as.character(response_cat),
    category == "Education" ~ as.character(response_cat),
    category == "Sexuality" ~ as.character(item),
    category == "Race/Ethnicity" ~ as.character(item),
    category == "Transgender & Gender-Diverse" ~ as.character(item),
    category == "Age" ~ as.character(item)
  )) |>
  select(-(starts_with("response")), -item) |>
  mutate(group = group |>
    replace_na("Prefer not to answer / Missing data") |>
    recode_factor(
      "Prefer not to answer" = "Prefer not to answer / Missing data",
      "prefer not to answer" = "Prefer not to answer / Missing data",
      "Missing Data" = "Prefer not to answer / Missing data"
    )
  ) |>
  summarise(.by = c(category, group), total = n_distinct(participant_id))

dem_totals <- d_demographics |>
  group_by() |>
  summarise(.by = category, total = sum(total)) |>
  mutate(group = "Total")

d_demographics <- d_demographics |>
  bind_rows(dem_totals) |>
  arrange(category, group)


# # Export----
write_csv(d_survey, "data/exp2_survey.csv")
write_csv(d_agg, "data/exp2_participant_covariates.csv")
write_csv(d_demographics, "data/exp2_demographics.csv")

# Loads accuracy data, sets up contrast coding and scaling----
exp3_load_data_acc <- function() {
  library(dplyr)
  library(forcats)
  library(scales)

  d <- read.csv("data/exp3_pronouns.csv", stringsAsFactors = TRUE) |>
    select(ParticipantID, Nametag, Intro, Pronoun_Pair, T_ID, Accuracy)

  # Remove trials with no pronouns
  d <- d |> filter(!is.na(Accuracy))

  # Mean-center effects code Nametag and Intro
  d$Nametag <- factor(d$Nametag, labels = c("-Nametag", "+Nametag"))
  contrasts(d$Nametag) <- cbind(c(-.5, .5))

  d$Intro <- factor(d$Intro, labels = c("-Intro", "+Intro"))
  contrasts(d$Intro) <- cbind(c(-.5, .5))

  # Orthogonal Helmert contrast codes for Pronoun Pair
  d <- d |> rename("Pronoun" = "Pronoun_Pair")
  d$Pronoun <- d$Pronoun |>
    fct_relevel("T_HS", after = 0) |>
    fct_relevel("HS_T", after = 1)
  contrasts(d$Pronoun) <- cbind(
    "Target" = c(-.66, +.33, +.33),
    "Dist"   = c(0,    -.50, +.50)
  )

  # Add dummy-coded factor for They vs He/She
  d <- d |> mutate(Pronoun_They0 = ifelse(Pronoun == "T_HS", 0, 1))

  # Dummy code Nametag and Intro
  d <- d |> mutate(
    Nametag_Yes0 = ifelse(Nametag == "+Nametag", 0, 1),
    Nametag_No0 = ifelse(Nametag == "-Nametag", 0, 1),
    Intro_Yes0 = ifelse(Intro == "+Intro", 0, 1),
    Intro_No0 = ifelse(Intro == "-Intro", 0, 1)
  )

  # Scale character (1-18)
  d <- d |> mutate(.keep = c("unused"), Character = rescale(T_ID, c(-0.5, 0.5)))

  # Subset and order
  d <- d |> select(
    ParticipantID, Nametag, Nametag_Yes0, Nametag_No0,
    Intro, Intro_Yes0, Intro_No0, Pronoun, Pronoun_They0, Character, Accuracy
  )

  return(d)
}

# Adds participant covariates to accuracy data, mean-centers + rescales them----
exp3_load_data_subj <- function() {
  # Join participant covariates to accuracy df
  d <- left_join(
      exp3_load_data_acc(),
      read.csv("data/exp3_participant-covariates.csv", stringsAsFactors = TRUE),
      by = "ParticipantID"
    ) |>
    rename("Familiarity" = "UseThey", "Rating" = "Rating_Name")

  # Remove participants with no pronouns (1) or no survey data (3)
  d <- d |> filter(!is.na(Age))
  d$ParticipantID <- droplevels(d$ParticipantID)

  # Scale THEN mean-center (on accuracy df)
  d <- d |> mutate(
    Age_C = scale(Age / 80, center = TRUE, scale = FALSE),
    Familiarity_C = scale(Familiarity / 2, center = TRUE, scale = FALSE),
    GenderBeliefs_C = scale(GenderBeliefs / 60, center = TRUE, scale = FALSE),
    LGBTQ_C = LGBQ - 0.50,
    Rating_C = scale(Rating / 6, center = TRUE, scale = FALSE),
    Sharing_C = scale(Sharing / 20, center = TRUE, scale = FALSE)
  )

  # Effects-code LGBTQ
  d <- d |> mutate(LGBTQ_Fct = as.factor(LGBTQ_C))
  contrasts(d$LGBTQ_Fct) <- cbind(c(-0.5, +0.5))

  # Subset and order
  d <- d |> select(
    ParticipantID, Nametag, Intro, Pronoun, Character, Accuracy,
    Age, Age_C, Familiarity, Familiarity_C,
    GenderBeliefs, GenderBeliefs_C, LGBTQ_C, LGBTQ_Fct,
    Rating, Rating_C, Sharing, Sharing_C
  )

  return(d)
}

# split into haves and create Pronoun effect vars for each half
exp3_d_reliability <- exp3_load_data_acc() |>
  select(ParticipantID, Pronoun, Accuracy) |>
  arrange(ParticipantID, Pronoun) |>  # sort by pronoun within participant
  mutate(Obs_Num = seq(1, length(Pronoun)))  |>
  mutate(Obs_Half = case_when(  # count odd and even trials
    is_even(Obs_Num) ~ "even",
    is_odd(Obs_Num)  ~ "odd"
  )) |>
  mutate(
    Pronoun_Even = case_when(  # effect of pronoun just in even trials
      Obs_Half == "even" & Pronoun == "T_HS" ~ -0.66,
      Obs_Half == "even" & Pronoun != "T_HS" ~ +0.33,
      Obs_Half == "odd" ~ 0
    ),
    Pronoun_Odd = case_when(  # effect of pronoun just in odd trials
      Obs_Half == "odd" & Pronoun == "T_HS" ~ -0.66,
      Obs_Half == "odd" & Pronoun != "T_HS" ~ +0.33,
      Obs_Half == "even" ~ 0
  ))

# run Bayesian model
exp3_m_reliability <- brm(
  formula = Accuracy ~ Pronoun_Even + Pronoun_Odd +  # fixed effects for halves
    (1 + Pronoun_Even + Pronoun_Odd | ParticipantID),  # random slopes by subj
  data = exp3_d_reliability,
  family = bernoulli(),  # keep default priors
  seed = 4, cores = 4,
  chains = 4, iter = 4000,
  file = "r_data/exp3_reliability"  # won't rerun because results are copied in
)
exp3_m_reliability

# tidy results
exp3_r_reliability <- exp3_m_reliability |>
  tidy() |>
  filter(str_detect(term, "Even") & str_detect(term, "Odd")) |>
  select(estimate, std.error, conf.low, conf.high) |>
  mutate(across(everything(), ~format(., digits = 2, nsmall = 2)))

exp3_r_reliability

# Run in parallel with 6 clusters
# Won't work when running as background job, but otherwise much faster
cl6 <- makeCluster(6)  # make 6 clusters, keep default type
clusterEvalQ(cl6, "buildmer")  # check all packages are loaded to each cluster
clusterExport(cl6, "exp3_d_subjCov")  # check data is loaded to each cluster

exp3_m_subj_cov <- buildmer(
  formula = Accuracy ~ Pronoun * Nametag * Intro *  # allow all interactions
            Age_C * Familiarity_C * GenderBeliefs_C * LGBTQ_Fct +
            Rating_C * Sharing_C +
            (1 | ParticipantID) + (1 | Character),
  data = exp3_d_subjCov,
  family = binomial,
  buildmerControl = list(
    direction = c("order", "backward"),  # max then backwards elim (default)
    cl = cl6,
    args = list(glmerControl(optimizer = "bobyqa")),  # nlminbwrap had huge SE
    # require Pronoun * Nametag * Intro and both random intercepts
    # aka keep hypothesis testing model
    include =
      "Pronoun * Nametag * Intro + (1 | ParticipantID) + (1 | Character)"
    )
  )
stopCluster(cl6)
remove(cl6)